#include <cctype>
#include <iostream>
#include <iomanip>
#include <string>
#include <sstream>

#include "TFile.h"
#include "TTree.h"
#include "TProfile2D.h"
#include "TH2F.h"
#include "TVirtualFFT.h"
#include "TMath.h"
#include "TFitResult.h"
#include "TGraph.h"
#include "TF1.h"

#define SAVE_WAVEFORMS 1

using namespace std;

void HammingWindow(float *window_constants, int N = 256);
void LowPassFilter(float *filter_constants);
float ButterworthLP(float f_in, float f_cutoff, int order = 4);

void CheckTime(int exp, int run, int module, 
	       string& output_file, int nentries) 
{
  char fiber_filename[1024];
  char camac_filename[1024];

  sprintf(fiber_filename,
	  "/home/craigb/belle/beamtest/pulser/exp99/e%02irun%04im%i.dst2.root",
	  exp,run,module);
  sprintf(camac_filename,
	  "/home/craigb/belle/beamtest/pulser/exp99/e%02irun%04i.camac",
	  exp,run);

  TFile *f1 = new TFile(fiber_filename);
  if (f1->IsZombie()) {
    cout << "Problem opening: " << fiber_filename << endl;
    return;
  }
  TTree *EventTree = (TTree *) f1->Get("EventTree");
  EventTree->AddFriend("camdata",camac_filename); // CAMAC is NOT your friend!!!

  short unsigned int ASIC_OriginWindow[4][4][8][4] = {{{{0}}}};
  float ASIC_ADC[4][4][8][4][64] = {{{{{0}}}}};
  float ftsw = 0; //This must be changed to an int for the original beamtest readout
  unsigned int qtm1 = 0;
  unsigned int qtm2 = 0;
  unsigned int EventNumber = 0;
  unsigned int camac_event = 0;
  unsigned short int raw_temperature = 0;
  EventTree->SetBranchAddress("ftsw",&ftsw);
  EventTree->SetBranchAddress("Temperature",&raw_temperature);
  //    EventTree->SetBranchAddress("qtm1",&qtm1);
  //    EventTree->SetBranchAddress("qtm2",&qtm2);
  EventTree->SetBranchAddress("ASIC_ADC",ASIC_ADC);
  EventTree->SetBranchAddress("ASIC_OriginWindow",ASIC_OriginWindow);
  EventTree->SetBranchAddress("EventNumber",&EventNumber);
  //    EventTree->SetBranchAddress("event",&camac_event);

  if (output_file.size() == 0) {
    stringstream ss; 
    ss << "e" << setfill('0') << setw(2) << exp;
    ss << "run" << setw(4) << run;
    ss << "m" << module << ".checktime_orig.root"; 
    output_file = ss.str();
  }

  cout << "Output will be written to `" << output_file << "'" << endl;

  TFile *temp_file = new TFile(output_file.c_str(),"RECREATE","stupid temporary file");
  char tree_name[256];
  sprintf(tree_name,"tree_m%i",module);
  TTree *tree = new TTree(tree_name,tree_name);
  float hit_time0[4][4][8] = {{{0}}};
  float hit_time1[4][4][8] = {{{0}}};
  float hit_time[4][4][8] = {{{0}}};
  float cfd_time0[4][4][8] = {{{0}}};
  float cfd_time1[4][4][8] = {{{0}}};
  float cfd_time[4][4][8] = {{{0}}};
  float ftsw_time = 0.0;
  float temperature = 0.0;
  float qtm_time = 0.0;
  float peak_adc[4][4][8] = {{{0}}};
  float premin_peak[4][4][8] = {0.0};
  float postmin_peak[4][4][8] = {{{0}}};
  float min_adc[4][4][8] = {{{0}}};
  int peak_sample[4][4][8] = {{{0}}};
  int min_sample[4][4][8] = {{{0}}};
  int threshold_sample0[4][4][8] = {{{0}}};
  int threshold_sample1[4][4][8] = {{{0}}};
  float threshold_sample[4][4][8] = {{{0}}};
  float threshold_slope[4][4][8] = {{{0}}};
  float cfd_sample[4][4][8] = {{{0}}};
  int channel[4][4][8] = {{{0}}};
  int hit_row[4][4][8] = {{{0}}};
  int hit_col[4][4][8] = {{{0}}};
  int quality[4][4][8] = {{{0}}};
  unsigned int LocalEventNumber = 0;
  float time[4][4][8][256] = {{{{0}}}};
  float adc[4][4][8][256] = {{{{0}}}};
  float power_above_250[4][4][8] = {{{0}}};
  float total_power[4][4][8] = {{{0}}};
  float fft_magnitude[4][4][8][129] = {{{{0}}}};
  float fft_cumul_sum[4][4][8][128] = {{{{0}}}};
  float fft_cumul_sum_fit_m[4][4][8] = {{{0}}};
  float fft_cumul_sum_fit_b[4][4][8] = {{{0}}};
  float fft_cumul_sum_fit_chi2[4][4][8] = {{{0}}};
  tree->Branch("temperature",&temperature,"temperature/F");
  tree->Branch("hit_time0",hit_time0,"hit_time0[4][4][8]/F");
  tree->Branch("hit_time1",hit_time1,"hit_time1[4][4][8]/F");
  tree->Branch("hit_time",hit_time,"hit_time[4][4][8]/F");
  tree->Branch("threshold_sample0",threshold_sample0,"threshold_sample0[4][4][8]/I");
  tree->Branch("threshold_sample1",threshold_sample1,"threshold_sample1[4][4][8]/I");
  tree->Branch("threshold_sample",threshold_sample,"threshold_sample[4][4][8]/F");
  tree->Branch("threshold_slope",threshold_slope,"threshold_slope[4][4][8]/F");
  tree->Branch("cfd_sample",cfd_sample,"cfd_sample[4][4][8]/F");
  tree->Branch("channel",channel,"channel[4][4][8]/i");
  tree->Branch("cfd_time0",cfd_time0,"cfd_time0[4][4][8]/F");
  tree->Branch("cfd_time1",cfd_time1,"cfd_time1[4][4][8]/F");
  tree->Branch("cfd_time",cfd_time,"cfd_time[4][4][8]/F");
  tree->Branch("peak_adc",peak_adc,"peak_adc[4][4][8]/F");
  tree->Branch("premin_peak",premin_peak,"premin_peak[4][4][8]/F");
  tree->Branch("postmin_peak",postmin_peak,"postmin_peak[4][4][8]/F");
  tree->Branch("min_adc",min_adc,"min_adc[4][4][8]/F");
  tree->Branch("peak_sample",peak_sample,"peak_sample[4][4][8]/I");
  tree->Branch("min_sample",min_sample,"min_sample[4][4][8]/I");
  tree->Branch("ftsw_time",&ftsw_time,"ftsw_time/F");
  tree->Branch("qtm_time",&qtm_time,"qtm_time/F");
  tree->Branch("LocalEventNumber",&LocalEventNumber,"LocalEventNumber/i");
  tree->Branch("EventNumber",&EventNumber,"EventNumber/i");
  //  tree->Branch("camac_event",&camac_event,"camac_event/i");
  tree->Branch("channel",channel,"channel[4][4][8]/I");
  tree->Branch("hit_row",hit_row,"hit_row[4][4][8]/I");
  tree->Branch("hit_col",hit_col,"hit_col[4][4][8]/I");
  tree->Branch("quality",quality,"quality[4][4][8]/I");
  tree->Branch("power_above_250",power_above_250,"power_above_250[4][4][8]/F");
  tree->Branch("total_power",total_power,"total_power[4][4][8]/F");
  tree->Branch("fft_cumul_sum_fit_m",fft_cumul_sum_fit_m,"fft_cumul_sum_fit_m[4][4][8]/F");
  tree->Branch("fft_cumul_sum_fit_b",fft_cumul_sum_fit_b,"fft_cumul_sum_fit_b[4][4][8]/F");
  tree->Branch("fft_cumul_sum_fit_chi2",fft_cumul_sum_fit_chi2,"fft_cumul_sum_fit_chi2[4][4][8]/F");
#if SAVE_WAVEFORMS
  tree->Branch("adc",adc,"adc[4][4][8][256]/F");
  tree->Branch("time",time,"time[4][4][8][256]/F");
  tree->Branch("fft_magnitude",fft_magnitude,"fft_magnitude[4][4][8][129]/F");
  tree->Branch("fft_cumul_sum",fft_cumul_sum,"fft_cumul_sum[4][4][8][128]/F");
#endif
  float threshold = -100;

  //Set up a TVirtualFFT for later
  int fft_size = 256;
  TVirtualFFT *fft  = TVirtualFFT::FFT(1,&fft_size,"R2C ES K");
  TVirtualFFT *ifft = TVirtualFFT::FFT(1,&fft_size,"C2R M K");

  float times[256];
  float dt_nominal  = 1./2.7;
  float f_refresh_clock = 127.216 / 6000.0;
  float T_refresh_clock = 1.0/f_refresh_clock;
  for (int i = 0; i < 128; ++i) {
    times[i] = dt_nominal * (float) i;
    times[i+128] = dt_nominal * (float) i + T_refresh_clock;
  }
  for (int col = 0; col < 4; ++col) {
    for (int row = 0; row < 4; ++row) {
      for (int ch = 0; ch < 8; ++ch) {
        channel[col][row][ch] = ch;
        hit_row[col][row][ch] = row;
        hit_col[col][row][ch] = col;
        quality[col][row][ch] = 0;
        for (int samp = 0; samp < 256; ++samp) {
          time[col][row][ch][samp] = times[samp];
          //This is some b.s. for reordering stuff
          //                                    if (samp % 2 == 0) {
          //                                            if (ch <= 5) {
          //                                                    time[col][row][ch][samp] += dt_nominal*2.0;
          //                                            } else {
          //                                                    time[col][row][ch][samp] += dt_nominal*4.0;
          //                                            }
          //                                    }
        }
      }
    }
  }

  //This array tracks the estimated order that samples should be placed into the array
  //It should correspond to the time ordering used above
  int placement[8][256];
  for (int ch = 0; ch < 8; ++ch) {
    for (int samp = 0; samp < 256; ++samp) {
      placement[ch][samp] = samp;
      if (samp % 2 == 0) {
        if (ch <= 3) {
          placement[ch][samp] -= 2;
        } else {
          placement[ch][samp] -= 4;
        }
      }
    }
  }

  //This array keeps track of the window function to be used for the FFT
  float window_function[256] = {0};
  HammingWindow(window_function);
  //Arrays to track the filter constants, if desired
  float re_fft_filter[256+1] = {0};
  LowPassFilter(re_fft_filter);
  //A fit function for linear fitting
  TF1 *linear = new TF1("linear","pol1",0,2.7);
  //A graph to perform fits to the cumulative FFT sum
  TGraph *fit_graph = new TGraph(42);

  if ((nentries == 0) || (nentries > EventTree->GetEntries()))
    nentries = EventTree->GetEntries();

  bool reported = false;
  for (int i = 0; i < nentries; ++i) {
    //TEMPORARILY KEEPING FILES SHORT FOR DEBUGGING
    //          if (i > 5000) {
    //                  break;
    //          }
    ////////////////////////
    LocalEventNumber = i;
    EventTree->GetEntry(i);
    if (EventNumber != camac_event) {
      if (!reported) {
        // cout << "Fiber event number was : " << EventNumber << endl;
        // cout << "CAMAC event number was : " << camac_event << endl;
        reported = true;
      }
    }

    //Rescale FTSW time to ns (from counts)
    ftsw_time = ftsw*0.025;
    qtm_time = (qtm1 + qtm2)/2.0 * 0.025;
    //          qtm_time = qtm1 * 0.025;
    //Rescale temperature to degrees C
    temperature = float(raw_temperature) * 0.0625;

    for (int col = 0; col < 4; ++col) {
      for (int row = 0; row < 4; ++row) {
        for (int ch = 0; ch < 8; ++ch) {
          quality[col][row][ch] = 0;
          float all_values[256] = {0};
          float all_values_reordered[256] = {0};
          for (int win = 0; win < 4; ++win) {
            for (int samp = 0; samp < 64; ++samp) {
              float adc_value;
              if (1) {
                //                                      if (ASIC_OriginWindow[col][row][ch][win] != 32) {
                adc_value = ASIC_ADC[col][row][ch][win][samp];
                //                                              adc_value -= ASIC_ADC[col][row][ch][win][(samp+1)%8];
              } else {
                adc_value = 0;
              }
              //                                        if (samp % 2 == 0) {
              //                                                adc_value *= 1.0;
              //                                        }
              all_values[win*64+samp] = adc_value;
              adc[col][row][ch][win*64+samp] = adc_value;
            }
          }
          //Reorder samples appropriately
          for (int j = 0; j < 256; ++j) {
            if (placement[ch][j] >= 0 && placement[ch][j] <= 255) {
              all_values_reordered[j] = all_values[ placement[ch][j] ];
            } else {
              all_values_reordered[j] = 0.0;
            }
          }
          for (int j = 0; j < 256; ++j) {
            all_values[j] = all_values_reordered[j];
            adc[col][row][ch][j] = all_values_reordered[j];
          }

          //////////////////
          //Apply the window function
          float windowed_points[256];
          for (int j = 0; j < 256; ++j) {
            windowed_points[j] = all_values_reordered[j] * window_function[j];
          }

          //Perform FFT and record power
          double fft_inout[2*((256+1)/2+1)];
          for (int point = 0; point < 256; ++point) {
            //                          fft_inout[point] = all_values_reordered[point];  //Uses the raw points (square window)
            fft_inout[point] = windowed_points[point];  //Uses the window function
          }
          fft->SetPoints(fft_inout);
          fft->Transform();
          double fft_re_part[256+1];
          double fft_im_part[256+1];
          fft->GetPointsComplex(fft_re_part,fft_im_part);
          float fft_mag[(256+1)];
          float fft_phi[(256+1)];
          for (int fft_points = 0; fft_points < 256+1; fft_points++) {
            fft_mag[fft_points] = sqrt( pow(fft_re_part[fft_points],2) +
                                        pow(fft_im_part[fft_points],2) );
            fft_phi[fft_points] = atan2( fft_im_part[fft_points] , fft_re_part[fft_points] );
            //Filter the magnitudes
            fft_mag[fft_points] *= re_fft_filter[fft_points];
            //Leave the phases untouched for now
          }
          //Only iterate over the half of the points we plan to draw into the tree
          for (int fft_points = 0; fft_points < 129; fft_points++) {
            fft_magnitude[col][row][ch][fft_points] = fft_mag[fft_points];
          }
          //Set up the fft_inout points for the inverse transform
          for (int fft_points = 0; fft_points < 256+1; fft_points++) {
            fft_re_part[fft_points] = fft_mag[fft_points] * cos(fft_phi[fft_points]);
            fft_im_part[fft_points] = fft_mag[fft_points] * sin(fft_phi[fft_points]);
          }
          ifft->SetPointsComplex(fft_re_part,fft_im_part);
          ifft->Transform();
          ifft->GetPoints(fft_inout);
          //Assign them to the original points to get the filtered waveform
          /*
            for (int j = 0; j < 256; ++j) {
            all_values_reordered[j] = fft_inout[j] / (256.+1.);
            //Unapply the window function
            all_values_reordered[j] /= window_function[j];
            adc[col][row][ch][j] = all_values_reordered[j];
            }
          */
          //Integrate power from 250 to 550 MHz
          float freq_bin_size = 2.7/2.0/128.0;
          float freq_start = 0.250;
          float freq_stop  = 0.550;
          float freq_bin_start = int(freq_start/freq_bin_size);
          float freq_bin_stop  = int(freq_stop/freq_bin_size);
          power_above_250[col][row][ch] = 0;
          total_power[col][row][ch] = 0;
          for (int bin = 0; bin < 128-1; ++bin) {
            float this_bin_integral = freq_bin_size * (fft_mag[bin]/2.0+fft_mag[bin+1]/2.0);
            if (bin > freq_bin_start && bin < freq_bin_stop) {
              power_above_250[col][row][ch] += this_bin_integral;
            }
            total_power[col][row][ch] += this_bin_integral;
            if (bin == 0) {
              fft_cumul_sum[col][row][ch][bin] = this_bin_integral;
            } else {
              fft_cumul_sum[col][row][ch][bin] = fft_cumul_sum[col][row][ch][bin-1] + this_bin_integral;
            }
          }
          fft_cumul_sum[col][row][ch][127] = fft_cumul_sum[col][row][ch][126];

          //Perform a fit to the cumulative fft power
          for (int bin = 0; bin < 42; ++bin) {
            int offset = 10;
            fit_graph->SetPoint(bin,freq_bin_size * float(bin+offset+1),fft_magnitude[col][row][ch][bin+offset]);
          }
          TFitResultPtr result = fit_graph->Fit("linear","QNS");
          fft_cumul_sum_fit_m[col][row][ch] = result->Value(1);
          fft_cumul_sum_fit_b[col][row][ch] = result->Value(0);
          fft_cumul_sum_fit_chi2[col][row][ch] = result->Chi2();

          //////////////////


          //Find min and max
          peak_adc[col][row][ch] = 0;
          min_adc[col][row][ch] = 0;
          peak_sample[col][row][ch] = -5000;
          min_sample[col][row][ch] = 5000;
          for (int j = 0; j < 256; ++j) {
            if (all_values[j] > peak_adc[col][row][ch]) {
              peak_adc[col][row][ch] = all_values[j];
              peak_sample[col][row][ch] = j;
            }
            if (all_values[j] < min_adc[col][row][ch]) {
              min_adc[col][row][ch] = all_values[j];
              min_sample[col][row][ch] = j;
            }
          }
          //Find the highest value in the N samples preceding the minimum
          int premin_scan_length = 10;
          int postmin_scan_length = 10;
          premin_peak[col][row][ch] = -5000;
          postmin_peak[col][row][ch] = -5000;
          for (int j = min_sample[col][row][ch]; j > min_sample[col][row][ch] - premin_scan_length && j >= 0; --j) {
            if ( all_values[j] > premin_peak[col][row][ch] ) {
              premin_peak[col][row][ch] = all_values[j];
            }
          }
          //Find the highest value in the N samples following the minimum
          for (int j = min_sample[col][row][ch]; j < min_sample[col][row][ch] + postmin_scan_length && j < 256; ++j) {
            if ( all_values[j] > postmin_peak[col][row][ch] ) {
              postmin_peak[col][row][ch] = all_values[j];
            }
          }
          //Find timing
          hit_time0[col][row][ch] = 0;
          hit_time1[col][row][ch] = 0;
          hit_time[col][row][ch] = 0;
          cfd_time0[col][row][ch] = 0;
          cfd_time1[col][row][ch] = 0;
          cfd_time[col][row][ch] = 0;
          //Simple threshold versions
          for (int j = min_sample[col][row][ch] - 20; j < 256; j += 2) { //Starts from 20 before min
            if (j % 2 != 0) j -= 1;
            if (j < 0) j = 0;
            //                  for (int j = 2; j < 256; j += 2) {  //Starts from the beginning
            if (all_values[j] < threshold) {
              //Simple version with no interpolation
              threshold_sample0[col][row][ch] = j-2;
              hit_time0[col][row][ch] = time[col][row][ch][j-2];
              //This part adds the interpolation
              float slope = (all_values[j] - all_values[j-2]) / (time[col][row][ch][j]-time[col][row][ch][j-2]);
              hit_time0[col][row][ch] += (threshold - all_values[j-2]) / slope;
              break;
            }
          }
          for (int j = min_sample[col][row][ch] - 20; j < 256; j += 2) { //Starts from 20 before min
            if (j % 2 == 0) j -= 1;
            if (j < 0) j = 1;
            //                  for (int j = 3; j < 256; j += 2) { //Starts from the beginning
            if (all_values[j] < threshold) {
              //Simple version with no interpolation
              threshold_sample1[col][row][ch] = j-2;
              hit_time1[col][row][ch] = time[col][row][ch][j-2];
              //This part adds the interpolation
              float slope = (all_values[j] - all_values[j-2]) / (time[col][row][ch][j]-time[col][row][ch][j-2]);
              hit_time1[col][row][ch] += (threshold - all_values[j-2]) / slope;
              break;
            }
          }
          for (int j = min_sample[col][row][ch] - 20; j < 256; ++j) { //Starts from 20 before min
            if (j < 0) j = 0;
            //                  for (int j = 1; j < 256; j += 1) {  //Starts from beginning
            if (all_values_reordered[j] < threshold) {
              //Simple version with no interpolation
              hit_time[col][row][ch] = time[col][row][ch][j-1];
              threshold_sample[col][row][ch] = j-1;
              //This part adds the interpolation
              float slope = (all_values_reordered[j] - all_values_reordered[j-1]) / (time[col][row][ch][j]-time[col][row][ch][j-1]);
              float slope_samples = (all_values_reordered[j] - all_values_reordered[j-1]);
              hit_time[col][row][ch] += (threshold - all_values_reordered[j-1]) / slope;
              threshold_sample[col][row][ch] += (threshold - all_values_reordered[j-1]) / slope_samples;

              threshold_slope[col][row][ch] = slope_samples;
              break;
            }
          }

	  /*
          //Constant fraction version
          int scan_point_odd  = (min_sample[col][row][ch] % 2 == 0) ? min_sample[col][row][ch]+1 : min_sample[col][row][ch] ;
          int scan_point_even = (min_sample[col][row][ch] % 2 == 0) ? min_sample[col][row][ch]   : min_sample[col][row][ch]+1 ;
          int scan_point = min_sample[col][row][ch];
          bool found_cf_odd = false;
          bool found_cf_even = false;
          bool found_cf = false;
          float fraction = 0.15;
          float target = fraction * min_adc[col][row][ch];
          while(!found_cf_odd) {
            if ( all_values[scan_point_odd-2] > target ) {
              float slope = (all_values[scan_point_odd] - all_values[scan_point_odd-2])/(time[col][row][ch][scan_point_odd] - time[col][row][ch][scan_point_odd-2]);
              float diff = target - all_values[scan_point_odd-2];
              cfd_time0[col][row][ch] = diff / slope + time[col][row][ch][scan_point_odd-2];
              found_cf_odd = true;
            }
            if (scan_point_odd > 3) {
              scan_point_odd -= 2;
            } else {
              break;
            }
          }
          while(!found_cf_even) {
            if ( all_values[scan_point_even-2] > target ) {
              float slope = (all_values[scan_point_even] - all_values[scan_point_even-2])/(time[col][row][ch][scan_point_even] - time[col][row][ch][scan_point_even-2]);
              float diff = target - all_values[scan_point_even-2];
              cfd_time1[col][row][ch] = diff / slope + time[col][row][ch][scan_point_even-2];
              found_cf_even = true;
            }
            if (scan_point_even > 2) {
              scan_point_even -= 2;
            } else {
              break;
            }
          }
          while(!found_cf) {
            if ( all_values_reordered[scan_point-1] > target ) {
              float slope = (all_values_reordered[scan_point] - all_values_reordered[scan_point-1])/(time[col][row][ch][scan_point] - time[col][row][ch][scan_point-1]);
              float slope_sample = (all_values_reordered[scan_point] - all_values_reordered[scan_point-1]);
              float diff = target - all_values_reordered[scan_point-1];
              cfd_time[col][row][ch] = diff / slope + time[col][row][ch][scan_point-1];
              cfd_sample[col][row][ch] = diff / slope_sample + scan_point;
              found_cf = true;
            }
            if (scan_point > 1) {
              scan_point -= 1;
            } else {
              break;
            }
          }
	  */
          //Implement "hard" cuts here
          //Junk data from FTSW.  Sometimes we get something well beyond
          //what a 12-bit number should look like.  We also lose events to
          //underflow when the trigger rate is too high.
          if ( ftsw_time < 10 || ftsw_time > 100 ) {
            quality[col][row][ch] |= 0x1;
          }
          //Hit times from waveforms are initialized to 0.  This should cut
          //out events where no hit past the given threshold was found.
          //We also filter out hits here that were outside the known possible time window.
          if ( min_adc[col][row][ch] > threshold || hit_time[col][row][ch] < 0. || hit_time[col][row][ch] > 110.) {
            quality[col][row][ch] |= 0x2;
          }
          //Window 32 has digital issues.  These manifest as very large
          //negative values for "min_adc".  These could be salvaged possibly
          //by masking out window 32 explicitly, but this cut for the moment
          //should remove events entirely where this occurred.
          if (min_adc[col][row][ch] < -1000) {
            quality[col][row][ch] |= 0x4;
          }
          //Maximum adc value should not be greater than the minimum one.  This
          //is indicative of at least one kind of cross talk.  This cut can be made
          //more severe as necessary.
          if (peak_adc[col][row][ch] + min_adc[col][row][ch] > 0) {
            quality[col][row][ch] |= 0x8;
          }
          //Peak after the minimum should be relatively large due to
          //overshoot/ringing effects.
          if (postmin_peak[col][row][ch] < 10 ) {
            quality[col][row][ch] |= 0x10;
          }
          //Peak before the minimum should be relatively small.  Cross talk hits
          //fluctuate high before they fluctuate low.
          if (premin_peak[col][row][ch] > 15 ) {
            quality[col][row][ch] |= 0x20;
          }
          //
        }
      }
    }
    temp_file->cd();
    tree->Fill();
  }
  tree->SetMarkerStyle(6);
  tree->AutoSave();
  temp_file->Write();

  { 
    float mypremin_peak[4][4][8];
    TFile* f = new TFile(output_file.c_str(), "READ");
    TTree* t = (TTree*) f->Get("tree_m0");
    
    t->SetBranchAddress("premin_peak", mypremin_peak);

    t->GetEntry(0);

    cout << "premin_peak[0][0][:] = [";
    for (int k=0; k < 8; k++) 
      cout << premin_peak[0][0][k] << ", ";
    cout << endl;

    cout << "min_adc[0][0][:] = [";
    for (int k=0; k < 8; k++) 
      cout << min_adc[0][0][k] << ", ";
    cout << endl;

    // cout << "mypremin_peak[0][0][:]  = [";
    // for (int k=0; k < 8; k++) 
    //   cout << mypremin_peak[0][0][k] << ", ";
    // cout << endl;
  }


}

void HammingWindow(float *window_constants, int N) {
  for (int i = 0; i < N; ++i) {
    window_constants[i] = 0.54 - 0.46 * cos(TMath::TwoPi() * float(i) / float(N-1));
  }
}

void LowPassFilter(float *re_filter_constants) {
  int N = 256+1;
  float freq_bin_size = 2.7/2.0/128.0;
  float freq_cutoff = 0.550;
  for (int i = 0; i < N/2; i++) {
    float freq = freq_bin_size * float(i+1);
    re_filter_constants[i] = ButterworthLP(freq , freq_cutoff);
  }
  for (int i = N; i >= N/2; i--) {
    float freq = freq_bin_size * float(N - i + 1);
    re_filter_constants[i] = ButterworthLP(freq , freq_cutoff);
  }

}

float ButterworthLP(float f_in, float f_cutoff, int order) {
  return ( 1. / ( 1.0 + pow(f_in/f_cutoff,order) ) );
}


void usage() {
  cout << "usage: btchecktime_orig"
       << " [opts] expno runno [modno]" << endl;
  cout << left << setw(20) << "\t-n N" << setw(50)
       << "Process only N events" << endl;
  cout << left << setw(20) << "\t-o output.root" << setw(50)
       << "Write output to output.root" << endl;
}


int main(int argc, char* argv[]) {
  int c;

  string output_file = "";
  int nentries = 100;

  while ((c = getopt(argc, argv, "o:n:h")) != -1) {
    switch(c)
      {
      case 'n':
        nentries = atoi(optarg);
        break;
      case 'o':
        output_file = optarg;
        break;
      case 'h':
        usage();
        return -2;
      case '?':
        if (optopt == 'o' || optopt == 'n') {
          cout << "Option -" << optopt << " requires an argument." << endl;
          cout << "See `btchecktime_orig -h' for details." << endl;
        }
        else if (isprint(optopt)) {
          cout << "Unknown option -" << optopt << endl;
          cout << "See `btchecktime_orig -h' for details." << endl;
        }
        else {
          cout << "Unrecognized option character " << optopt << endl;
          cout << "See `btchecktime_orig -h' for details." << endl;
        }
        return 1;
      default:
        abort();
      }
  }

  if ((argc - optind) < 2) {
    cout << "Please supply an experiment number and a run number." << endl;
    cout << "See `btchecktime_orig -h' for details." << endl;
    return -2;
  }
 
  int expno = atoi(argv[optind++]);
  int runno = atoi(argv[optind++]);

  int modno = 0;  
  if (optind < argc)
    modno = atoi(argv[optind]);

  CheckTime(expno, runno, modno, output_file, nentries);

  return 0;
}

